using System;

namespace PVTlibrary
{
    [Serializable]
    public class PROPS
    {

        public PROPS()
        {

        }

        /// <summary>
        /// liquid heat capacity by Rowlinson/Bondi correlation
        /// </summary>
        /// <param name="cpig">ideal gas heat capacity, kJ/kg.K</param>
        /// <param name="T">temperature in K</param>
        /// <param name="Tc">critical temperature in K</param>
        /// <param name="w">acentric factor</param>
        /// <param name="MW">mole weight, kg/kmol</param>
        /// <returns></returns>
        double Cpl_rb(double cpig , double T , double Tc , double w , double MW ) 
        {

            var R = 8.314;
            var Tr = T / Tc;

            var cpl = cpig * MW + 1.45 * R + 0.45 * R * 1/(1 - Tr) + 0.25 * w * R * (17.11 + 25.2 * Math.Pow(1 - Tr, 1.0 / 3.0) * 1/Tr + 1.742 * 1/(1 - Tr) );

            return cpl / MW;// 'kJ/kg.K*/
        }


        /// <summary>
        /// Ideal gas state heat capacity for petroleum fractions (pseudos) by Lee-Kesler
        /// </summary>
        /// <param name="WK"></param>
        /// <returns></returns>
        double Cpig_lk(double WK, double w , double  T ) 
        {

            var A0 = -1.41779 + 0.11828 * WK;
            var A1 = -(6.99724 - 8.69326 * WK + 0.27715 * WK * WK) * 0.0001;
            var A2 = -0.0000022582;
            var B0 = 1.09223 - 2.48245 * w;
            var B1 = -(3.434 - 7.14 * w) * 0.001;
            var B2 = -(7.2661 - 9.2561 * w) * 0.0000001;
            var C = Math.Pow(((12.8 - WK) * (10 - WK) / (10 * w)),2);

            var Cp = A0 + A1 * T + A2 * T * T - C * (B0 + B1 * T + B2 * T *T) ;//'kJ/kg

            return Cp;
         }


        /// <summary>
        /// Estimativa da tensão interfacial pelo método de Brock e Bird
        /// </summary>
        /// <param name="T"></param>
        /// <returns></returns>
        double  sigma_bb(double T , double Tb , double Tc, double Pc ) 
        {

            //Dim Tr, Tbr, Q, tmp, Pc_ As Double

            var Pc_ = Pc / 100000 ;

            var Tr = T / Tc;

            var Tbr = Tb / Tc;

            var Q = 0.1196 * (1 + Tbr * Math.Log(Pc_ / 1.01325) / (1 - Tbr)) - 0.279;

            var tmp = Math.Pow(Pc_ , 2.0 / 3.0) * Math.Pow(Tc ,1.0 / 3.0) * Q * Math.Pow(1 - Tr,11.0 / 9.0);

            tmp = tmp / 1000; //'N/m

           return tmp;
        }

        /// <summary>
        /// 'Estimativa da viscosidade de líquidos pelo método de Letsou e Stiel
        /// </summary>
        /// <param name="T"></param>
        /// <returns></returns>
        double  viscl_letsti(double T, double Tc, double Pc, double w, double MM) 
        {           
            var Tr = T / Tc;

            Pc = Pc / 100000;

            var e0 = (2.648 - 3.725 * Tr + 1.309 * Tr * Tr) * 0.001;

            var e1 = (7.425 - 13.39 * Tr + 5.933 * Tr * Tr) * 0.001;

            var e = 0.176 * Math.Pow(Tc / (Math.Pow(MM , 3) * Math.Pow(Pc ,4)), 1.0 / 6.0);

            var viscl_letsti = (e0 + e1) / e / 1000 ;//'Pa.s

            return viscl_letsti;
        }

        /// <summary>
        /// Estimativa da viscosidade de gases pelo método de Lucas
        /// </summary>
        /// <param name="T"></param>
        /// <param name="Tc"></param>
        /// <param name="Pc"></param>
        /// <param name="w"></param>
        /// <param name="MM"></param>
        /// <returns></returns>
        double  viscg_lucas(double T, double Tc, double Pc, double w, double MM) 
        {
            var Tr = T / Tc;

            Pc = Pc / 100000;

            var e = 0.176 * Math.Pow(Tc / (Math.Pow(MM , 3) * Math.Pow(Pc,4)), 1.0 / 6.0);

            var  tmp = (1 / e) * (0.807 * Math.Pow(Tr, 0.618) - 0.357 * Math.Exp(-0.449 * Tr) + 0.34 * Math.Exp(-4.058 * Tr) + 0.018);

            var viscg_lucas = tmp / 1000000.0 * 100 / 1000 ;//'Pa.s

            return viscg_lucas;
        }


        /// <summary>
        /// Cálculo da densidade do líquido utilizando o método de Rackett
        /// 'Modified HBT method - Thomson (para líquidos comprimidos)
        /// </summary>
        /// <param name="T"></param>
        /// <param name="Tc"></param>
        /// <param name="Pc"></param>
        /// <param name="w"></param>
        /// <param name="MM"></param>
        /// <param name="ByVal"></param>
        /// <returns></returns>
        double liq_dens_rackett(double T, double Tc, double Pc, double w, double MM, double ZRa  = 0, double P = 0, double Pvp  = 0) 
        {          

            var liq_dens_rackett=0.0;
            if( double.IsNaN(Pvp))
                Pvp = 0;

            var R = 83.14;

            Pc = Pc / 100000;

            var Tr = T / Tc;

            if (ZRa == 0.0)
                ZRa = 0.29056 - 0.08775 * w;

            var tmp = R * Tc / Pc * Math.Pow(ZRa , 1 + Math.Pow(1 - Tr,2.0 / 7.0));

            if(Pvp != 0.0)
            {                
                var a = -9.070217;
                var b = 62.45326;
                var d = -135.1102;
                var f = 4.79594;
                var g = 0.250047;
                var h = 1.14188;
                var j = 0.0861488;
                var k = 0.0344483;

                var e = Math.Exp(f + g * w + h * w *w);
                var c = j + k * w;

                var  beta = Pc * 100000 * (-1 + a * (1 - Tr) ^ (1 / 3) + b * (1 - Tr) ^ (2 / 3) + d * (1 - Tr) + e * (1 - Tr) ^ (4 / 3));


                var tmpcor = tmp * (1 - c * Math.Log((beta + P) / (beta + Pvp)));

                liq_dens_rackett = 0.001 * MM / (tmpcor * 0.000001);// 'kg/m3 ''''m3/mol

            }
            else
                liq_dens_rackett = 0.001 * MM / (tmp * 0.000001);// 'kg/m3 ''''m3/mol

            return liq_dens_rackett;
        }


        /// <summary>
        /// 
        /// </summary>
        /// <param name="T"></param>
        /// <param name="Tpc"></param>
        /// <param name="Ppc"></param>
        /// <param name="Cp"></param>
        /// <param name="fluido"></param>
        /// <param name="SG"></param>
        /// <returns></returns>
        double JT_Goldzberg(double T, double Tpc, double Ppc , double Cp , string fluido , double SG )
        {
            //'T em °F
            T = 5 / 9 * (T + 459.7);
            //'Tpc em °F
            Tpc = 5 / 9 * (Tpc + 459.7);
            //'Ppc em lbf/in.2
            Ppc = Ppc / 6894.76;
            //'Cp em J/kg°F
            Cp = Cp * 1.8 / 1000;
            //'tipo como "V" ou "L"
            //'SG = densidade absoluta ou densidade do líquido em lbm/ft3

            //'JT em °F in.2/lbf

            //'para gases

            var Tpr = T / Tpc;
            var eta=0.0;

            if (fluido == "V") 
                eta = 0.0048823 * Tpc * (18 / Tpr ^ 2 - 1) / (Ppc * Cp * SG);

            //'para líquidos
            if (fluido == "L")  
                eta = -1 / (SG / 16.0185 * Cp);

            var JT_Goldzberg = eta / 1.8 / 6894.76;  //'(K/Pa)
            
            return JT_Goldzberg;
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="T"></param>
        /// <param name="Tb"></param>
        /// <param name="Tc"></param>
        /// <param name="M"></param>
        /// <param name="Tipo"></param>
        /// <returns></returns>
        double condl_latini(double T, double Tb, double Tc, double M, string Tipo )
        {
           // Dim Tr, A, A_, alpha, beta, lambda, gamma

            var Tr = T / Tc;

            //'Hidrocarbonetos saturados (Tipo = "H")

            var A_ = 0.0035;
            var alpha = 1.2;
            var beta = 0.5;
            var gamma = 0.167;

            if (Tipo == "O")
            {
                //'Oleofinas
                A_ = 0.0361;
                alpha = 1.2;
                beta = 1;
                gamma = 0.167;
            }
            else if(Tipo == "C") 
            {
                //'Cicloparafinas
                A_ = 0.031;
                alpha = 1.2;
                beta = 1.0;
                gamma = 0.167;
            }
            else if(Tipo == "A")
            {
                //'Aromáticos
                A_ = 0.0346;
                alpha = 1.2;
                beta = 1.0;
                gamma = 0.167;
            }
            else if(Tipo == "X" )
            {
                //'Outros (ex. água)
                A_ = 0.494;
                alpha = 0;
                beta = 0.5;
                gamma = -0.167;
            }

            var A = A_ * Math.Pow(Tb ,alpha) / (Math.Pow(M,beta) * Math.Pow(Tc,gamma));
            var lambda = A * (1 - Tr) ^ 0.38 / Tr ^ (1 / 6);

            if( T / Tc > 0.98)
                lambda = 0;

            var condl_latini = lambda ;  //'W/(m.K)*/

            return condl_latini;
        }

        double condtg_elyhanley(double T, double Tc, double Vc, double Zc, double w, double M, double Cv)
        {
           
            //Dim Tr, Tplus, teta, omega, f, h, T0, eta0, lambda0, C(8), lambda_, eta_, Hgrande
            double Tplus;
            double[] C = new double[9];

            var Tr = T / Tc;

            M = M / 1000;

            if(Tr <= 2)
                Tplus = Tr;
            else
                Tplus = 2;


            var teta = 1 + (w - 0.011) * (0.56553 - 0.86276 * Math.Log(Tplus) - 0.69852 / Tplus);

            var omega = (1 + (w - 0.011) * (0.3856 - 1.1617 * Math.Log(Tplus))) * 0.288 / Zc;

            var f = Tc * teta / 190.4;

            var h = Vc * omega / 99.2 * 1000000.0;

            var T0 = T / f;

            C[0] =  2907740.0;
            C[1] =  -3312870.0;
            C[2] =  1608100.0;
            C[3] =  -433190.0;
            C[4] =  70624.8;
            C[5] =  -7116.62;
            C[6] =  432.517;
            C[7] =  -14.4591;
            C[8] =  0.203712;

            int i = 0;

            var eta0 = 0.0;

            do 
            {
                eta0 = eta0 + 0.0000001 * C(i) * T0 ^ ((i - 3) / 3);
                i = i + 1;
            }
            while( i < 9);

            //'eta0 = viscg_lucas(T, 190.6, 45.4 * 101325, 0.012, 16)

            var lambda0 = 1944 * eta0;

            var Hgrande = (0.01604 / M) ^ 0.5 * f ^ 0.5 / h ^ (2 / 3);

            var lambda_ = lambda0 * Hgrande;

            var eta_ = eta0 * Hgrande * M / (0.01604);

            var condtg_elyhanley = lambda_ + (eta_ / M) * 1.32 * (Cv - 3 * 8.314 / 2) ;//'W/(m.K)
            

            return condtg_elyhanley;
        }

        double Tcm(double[]  Vz, double[]  VTc)
        {
            /*
            Dim i, n

            n = UBound(Vz)

            i = 0

            Tcm = 0

            Do

                Tcm += Vz(i) * VTc(i)

                i = i + 1

            Loop Until i = n + 1*/

            double Tcm=0.0;
            var n = Vz.Length;

            for(int i=0;i<n;i++)
            {
                Tcm += Vz[i] * VTc[i];
            }

            return Tcm;
        }

        double Pcm(double[]  Vz, double[] VPc)
        {

            double Pcm=0.0;
            var n = Vz.Length;

            for(int i=0;i<n;i++)
            {
                Pcm += Vz[i] * VPc[i];
            }

            return Pcm;
        }

        double Vcm(double[]  Vz, double[] VVc) 
        {
            double Vcm=0.0;
            var n = Vz.Length;

            for(int i=0;i<n;i++)
            {
                for(int j=0;j<n;j++)
                {
                    Vcm += Vz[i] * Vz[j] * (0.5 * VVc[i] ^ (1 / 3) + 0.5 * VVc[j] ^ (1 / 3)) ^ 3;
                }
            }

            return Vcm;
        }

        double Zcm(double[]  Vz, double[] VZc)
        {
            double Zcm=0.0;
            var n = Vz.Length;

            for(int i=0;i<n;i++)
            {
                Zcm += Vz[i] * VZc[i];
            }

            return Zcm;
        }

        double wm(double[]  Vz, double[] Vw) 
        {

            double wm=0.0;
            var n = Vz.Length;

            for(int i=0;i<n;i++)
            {
                wm += Vz[i] * Vw[i];
            }

            return wm;
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="Vz"></param>
        /// <param name="VMM"></param>
        /// <returns></returns>
        double MMm(double[]  Vz, double[]  VMM) 
        {
            double MMm=0.0;
            var n = Vz.Length;

            for(int i=0;i<n;i++)
            {
                MMm += Vz[i] * VMM[i];
            }

            return MMm;
        }

        double viscg_jossi_stiel_thodos(double eta0 , double T, double V, double Tc, double Pc, double Vc, double MM) 
        {
            //'Estimativa da viscosidade de gases em pressões elevadas pelo 
            //'método de Jossi-Stiel-Thodos

            var Tr = T / Tc;

            Pc = Pc / 101325;

            var rho_r = Vc / V;

            var xit = Math.Pow(Tc / Math.Pow(MM,3) / Math.Pow(Pc,4), 1.0 / 6.0);

            eta0 = eta0 * 1000000.0 / 100 * 1000;

            var tmp = ((1.023 + 0.23364 * rho_r + 0.58533 * Math.Pow(rho_r,2) - 0.40758 * Math.Pow(rho_r,3) + 0.093324 * Math.Pow(rho_r ,4) )  - 1) / xit + eta0;

            var viscg_jossi_stiel_thodos = tmp / 10000 / 1000 ;//'Pa.s

            return viscg_jossi_stiel_thodos;
        }

        ////
        double  Vc(double Tc, double Pc, double w)
        {
            var Vc = 8.314 * (0.291 - 0.08 * w) * Tc / Pc * 1000 ; //'m3/kmol
            return Vc;
        }

        double  Vc(double Tc, double Pc, double w, double Zc)
        {
            var Vc = 8.314 * Zc * Tc / Pc * 1000;  //'m3/kmol
            return Vc;
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="w"></param>
        /// <returns></returns>
        double Zc1(double w)
        {
            var Zc1 = 0.291 - 0.08 * w; 
            return Zc1;
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="VVl"></param>
        /// <param name="Vcondl"></param>
        /// <param name="Vz"></param>
        /// <returns></returns>
        double condlm_li(double[] VVl, double[]  Vcondl, double[]  Vz)
        {

            double somaz=0.0;
            var n = Vz.Length;

            for(int i=0;i<n;i++)
            {
                if(  Vz[i] != 0 &&  !double.IsNaN(VVl[i]))
                    somaz = somaz +  Vz[i] * VVl[i];
            }

            var phi = new double[n];

           for(int i=0;i<n;i++)
            {
                if(  Vz[i] != 0 &&  !double.IsNaN(VVl[i]) )
                    phi[i] =   Vz[i] * VVl[i] / somaz;
            }

           double condlm = 0;
           for(int i=0;i<n;i++)
            {
                for(int j=0;j<n;j++)
                {
                    if( Vz[i]!= 0.0 && Vz[j] != 0.0)  
                        condlm += phi[i] * phi[j] * (2 / (1.0/Vcondl[i] + 1.0/Vcondl[j]) );
                }
            }

            return condlm;    //'W/(m.K)
            
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="VCOEFF"></param>
        /// <param name="T"></param>
        /// <returns></returns>
        double CPig_pol(double[] VCOEFF, double T)
        {
            var T2  =    T*T;
            var T3  =    T2*T;
            var T4  =    T3*T;
            var tmp2 = VCOEFF[0] + VCOEFF[1] * T + VCOEFF[2] * T2 + VCOEFF[3] *T3 + VCOEFF[4] * T4;

            return tmp2;
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="T"></param>
        /// <param name="Tc"></param>
        /// <param name="Tb"></param>
        /// <param name="Pc"></param>
        /// <returns></returns>
        double Vc(double T, double Tc, double Tb, double Pc)
        {
            Pc = Pc / 101325;
            var R = 82.0575;
            var s = Tb / Tc * Math.Log(Pc) / (1 - Tb / Tc);
            var tmp = (R * Tc / Pc) * (0.3883 - 0.0179 * s);

            var Vc = tmp * 0.000001;
            return Vc;
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="T"></param>
        /// <param name="T1"></param>
        /// <param name="T2"></param>
        /// <param name="v1"></param>
        /// <param name="v2"></param>
        /// <returns></returns>
        double oilvisc_twu(double T, double T1, double T2, double v1, double v2)
        {
            //'v = m2/s, T = K

            //var Z, Z1, Z2, B, vk1, vk2;

            var vk1 = v1 * 1000000.0;
            var vk2 = v2 * 1000000.0;

            T = 1.8 * T;
            T1 = 1.8 * T1;
            T2 = 1.8 * T2;

            var  Z1 = vk1 + 0.7 + Math.Exp(-1.47 - 1.84 * vk1 - 0.51 * vk1 * vk1);
            var  Z2 = vk2 + 0.7 + Math.Exp(-1.47 - 1.84 * vk2 - 0.51 * vk2 * vk2);

            var var1 = (Math.Log(Math.Log(Z1)) + Math.Log(Math.Log(Z2)));
            var var2 = (Math.Log(T1) + Math.Log(T2));
            var B = var1 / var2;

            var Z = Math.Exp(Math.Exp(Math.Log(Math.Log(Z1)) + B * (Math.Log(T) - Math.Log(T2))));

            var  tmp = Z - 0.7 - Math.Exp(-0.7487 - 3.295 * (Z - 0.7) + 0.6119 * (Z - 0.7) * (Z - 0.7) - 0.3193 * (Z - 0.7) * (Z - 0.7)*(Z - 0.7));

            var oilvisc_twu = tmp * 0.000001;  //'m2/s
        
            return oilvisc_twu;
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="TIPO"></param>
        /// <param name="T"></param>
        /// <param name="P"></param>
        /// <param name="Vz"></param>
        /// <param name="VKij"></param>
        /// <param name="Vzmass"></param>
        /// <param name="VTc"></param>
        /// <param name="VPc"></param>
        /// <param name="VCpig"></param>
        /// <param name="VMM"></param>
        /// <param name="Vw"></param>
        /// <param name="VZRa"></param>
        /// <returns></returns>
        double CpCvR(string TIPO, double T, double P, double[] Vz, double[,] VKij, double Vzmass, double[] VTc, double[] VPc, double[] VCpig, double[] VMM, double[] Vw, double VZRa)
        {
                /*
            Dim ai(), bi(), ci() As Double
            Dim n, R As Double
            Dim Tc(), Pc(), Vc(), w(), Zc(), alpha(), m(), a(,), b(,), Z, Tr() As Double
            Dim i, j, dadT
            */
            var n = Vz.Length;

            var ai =  new double[n];
            var bi =  new double[n];
            var ci =  new double[n];
            var a = new double[n,n];
            var b = new double[n,n];
            var Tc =  new double[n];
            var Pc =  new double[n];
            var Vc =  new double[n];
            var Zc =  new double[n];
            var w =  new double[n];
            var alpha =  new double[n];
            var m =  new double[n];
            var Tr =  new double[n];

            var coeff = new double[4];

            var R = 8.314;

            for(int i=0;i<n;i++)
            {
                Tc[i] = VTc[i];
                Tr[i] = T / Tc[i];
                Pc[i] = VPc[i];
                w[i]  = Vw[i];
            }

            var MMm = 0;
            for(int i=0;i<n;i++)
            {
                MMm += Vz[i] * VMM[i];
            }


            for(int i=0;i<n;i++)
            {
                alpha[i] = (1 + (0.37464 + 1.54226 * w[i] - 0.26992 * w[i] ^ 2) * (1 - (T / Tc[i]) ^ 0.5)) ^ 2;
                ai[i] = 0.45724 * alpha[i] * R ^ 2 * Tc[i] ^ 2 / Pc[i];
                bi[i] = 0.0778 * R * Tc[i] / Pc[i];
                ci[i] = 0.37464 + 1.54226 * w[i] - 0.26992 * w[i] ^ 2;
            }


            for(int i=0;i<n;i++)
            {
                for(int j=0;j<n;j++)
                {
                    a[i, j] = (ai[i] * ai[j]) ^ 0.5 * (1 - VKij[i, j]);
                }
            }

            var am = 0.0;
            for(int i=0;i<n;i++)
            {
                for(int j=0;j<n;j++)
                {
                    am = am + Vz[i] * Vz[j] * a[i, j];
                }
            }

            var bm = 0.0;
            for(int i=0;i<n;i++)
            {
                bm = bm + Vz[i] * bi[i];
            }



            var AG1 = am * P / Math.Pow(R * T, 2);
            var BG1 = bm * P / (R * T);


            coeff[0] = -AG1 * BG1 + BG1 * BG1 + BG1 *BG1*BG1;
            coeff[1] = AG1 - 3 * BG1 ^ 2 - 2 * BG1;
            coeff[2] = BG1 - 1;
            coeff[3] = 1;

            var temp1 = Poly_Roots(coeff);
            double tv;
            double tv2;
            /*
            If Not IsNumeric(temp1) Then

                If temp1(0, 0) > temp1(1, 0) Then
                    tv = temp1(1, 0)
                    tv2 = temp1(1, 1)
                    temp1(1, 0) = temp1(0, 0)
                    temp1(0, 0) = tv
                    temp1(1, 1) = temp1(0, 1)
                    temp1(0, 1) = tv2
                End If
                If temp1(0, 0) > temp1(2, 0) Then
                    tv = temp1(2, 0)
                    temp1(2, 0) = temp1(0, 0)
                    temp1(0, 0) = tv
                    tv2 = temp1(2, 1)
                    temp1(2, 1) = temp1(0, 1)
                    temp1(0, 1) = tv2
                End If
                If temp1(1, 0) > temp1(2, 0) Then
                    tv = temp1(2, 0)
                    temp1(2, 0) = temp1(1, 0)
                    temp1(1, 0) = tv
                    tv2 = temp1(2, 1)
                    temp1(2, 1) = temp1(1, 1)
                    temp1(1, 1) = tv2
                End If

                If TIPO = "L" Then
                    Z = temp1(0, 0)
                    If temp1(0, 1) <> 0 Then
                        Z = temp1(1, 0)
                        If temp1(1, 1) <> 0 Then
                            Z = temp1(2, 0)
                        End If
                    End If
                    If Z < 0 Then Z = temp1(1, 0)
                ElseIf TIPO = "V" Then
                    Z = temp1(2, 0)
                    If temp1(2, 1) <> 0 Then
                        Z = temp1(1, 0)
                        If temp1(1, 1) <> 0 Then
                            Z = temp1(0, 0)
                        End If
                    End If
                End If

            Else

                Dim findZV, dfdz, zant As Double
                If TIPO = "V" Then Z = 1 Else Z = 0.05
                Do
                    findZV = coeff(3) * Z ^ 3 + coeff(2) * Z ^ 2 + coeff(1) * Z + coeff(0)
                    dfdz = 3 * coeff(3) * Z ^ 2 + 2 * coeff(2) * Z + coeff(1)
                    zant = Z
                    Z = Z - findZV / dfdz
                    If Z < 0 Then Z = 1
                Loop Until Math.Abs(findZV) < 0.0001 Or Double.IsNaN(Z)


            End If

            Dim V = (Z * R * T / P) ' m3/mol

            Dim tmp1 = MMm / V / 1000

            Dim aux1 = -R / 2 * (0.45724 / T) ^ 0.5
            i = 0
            Dim aux2 = 0
            Do
                j = 0
                Do
                    aux2 += Vz(i) * Vz(j) * (1 - VKij(i, j)) * (ci(j) * (ai(i) * Tc(j) / Pc(j)) ^ 0.5 + ci(i) * (ai(j) * Tc(i) / Pc(i)) ^ 0.5)
                    j = j + 1
                Loop Until j = n + 1
                i = i + 1
            Loop Until i = n + 1

            dadT = aux1 * aux2
            Dim d2adt2 = R / 4 * (0.45724 / T) ^ 0.5 * (1 / T) * aux2
            'Dim d2adt2 = 0.169049 * R / (T ^ (3 / 2))

            Dim dP_dT_V = R / (V - bm) - dadT / (V ^ 2 + 2 * bm * V - bm ^ 2)

            Dim dV_dT_P = dP_dT_V / (R * T / (V - bm) ^ 2 - am * (2 * V + 2 * bm) / (V * (V + bm) + bm * (V - bm)) ^ 2)

            Dim dP_dV_T = -R * T * (V - bm) ^ -2 - am * (V ^ 2 + 2 * bm * V - bm ^ 2) ^ -2 * (2 * V + 2 * bm)

            Dim d2P_dT2 = -1 / (V ^ 2 + 2 * bm * V - bm ^ 2) * d2adt2

            Dim var = (bm + V) / (2 ^ 0.5 * bm)

            Dim Int_d2P_dT2_V_dV = -d2adt2 * Math.Log((-(2 ^ 0.5) * bm + bm + V) / ((2 ^ 0.5) * bm + bm + V)) / (8 ^ 0.5 * bm)

            Dim Cpm_ig = 0
            i = 0
            Do
                Cpm_ig += Vzmass(i) * VCpig(i) * MMm
                i += 1
            Loop Until i = n + 1

            Dim Cv = T * Int_d2P_dT2_V_dV + Cpm_ig - 2 * R - T * dP_dT_V ^ 2 / dP_dV_T
            'Dim Cp = Cpm_ig + T * Int_d2P_dT2_V_dV - T * dP_dT_V ^ 2 / dP_dV_T - R
            Dim Cp = Cpm_ig - R + T * dP_dT_V * dV_dT_P - T * d2adt2 / (8 ^ 0.5 * bm) * Math.Log((V + (1 - 2 ^ 0.5) * bm) / (V + (1 + 2 ^ 0.5) * bm))

            Dim Cp_Cv2 = Cp / Cv

            Dim Cp_Cv = 1 - (T * dP_dT_V ^ 2 / dP_dV_T) / (Cpm_ig - R + T * Int_d2P_dT2_V_dV)
            'Cv = Cp / Cp_Cv

            Dim tmp(2) As Double
            tmp(0) = Cp_Cv2
            tmp(1) = Cp / MMm
            tmp(2) = Cv / MMm

            CpCvR = tmp*/
        }

        double  JT_PR(double TIPO, double T, double P, double[] Vz, double[] Vzmass, double VTc, double VPc, double VCpig, double VMM, double Vw, double VZRa)
        {
            /*
            //Dim n, R, Cpm_ig As Double
            var  vetor = new double[8];
            //Dim Tc(), Pc(), Vc(), W(), Zc(), a, b, c, Tr() As Double

            var n = Vz.length;

            //ReDim Zc(n), Tc(n), Pc(n), Vc(n), W(n), Tr(n)
            Zc  =new double[n];
            Tc  =new double[n];
            Pc  =new double[n];
            Vc  =new double[n];
            W   =new double[n];
            Tr  =new double[n];

            R = 8.314;

            Dim i, j As Integer
            i = 0
            Do
                Tc(i) = VTc(i)
                Tr(i) = T / Tc(i)
                Pc(i) = VPc(i)
                W(i) = Vw(i)
                Vc(i) = Zc1(W(i)) * R * Tc(i) / Pc(i)
                i = i + 1
            Loop Until i = n + 1

            i = 0
            Dim Vcm = 0
            Dim wm = 0
            Dim Zcm = 0
            Dim MMm = 0
            Dim ZRam = 0
            Do
                If Vz(i) <> 0 Then
                    Vcm += Vz(i) * Vc(i)
                    wm += Vz(i) * W(i)
                    Zcm += Vz(i) * Zc1(W(i))
                    MMm += Vz(i) * VMM(i)
                    ZRam = Vz(i) * VZRa(i)
                End If
                i += 1
            Loop Until i = n + 1

            i = 0
            Dim Tcm = 0
            Do
                j = 0
                Do
                    If Vz(i) <> 0 And Vz(j) <> 0 Then Tcm += Vz(i) * Vz(j) * (Tc(i) * Tc(j)) ^ 0.5
                    j += 1
                Loop Until j = n + 1
                i += 1
            Loop Until i = n + 1

            Dim Pcm = Zcm * R * Tcm / (Vcm)

            Dim V = 0
            If TIPO = "L" Then

                'V = (Z_PR(T, P, Vz, VTc, VPc, Vw, "L") * R * T / P) * 1000 ' m3/kgmol

            ElseIf TIPO = "V" Then

                'V = (Z_PR(T, P, Vz, VTc, VPc, Vw, "V") * R * T / P) * 1000 ' m3/kgmol

            End If

            a = 0.45724 * R ^ 2 * Tcm ^ 2 / Pcm
            b = 0.0778 * R * Tcm / Pcm
            c = 0.37464 + 1.54226 * wm - 0.26992 * wm ^ 2

            Dim Trm = T / Tcm
            Dim AG = -b * R * T
            Dim BG = -2 * a * T * (1 + c - c * Trm ^ 0.5) * (-0.5 * c * Trm ^ 0.5)
            Dim CG = a * (1 + c - c * Trm ^ 0.5) ^ 2

            Dim dP_dT_V = R / (V - b) - (a * (1 + c) * c * Tcm ^ -0.5 * T ^ -0.5 + a * c ^ 2 * Tcm ^ -1) / (V ^ 2 + 2 * b * V - b ^ 2)

            Dim dP_dV_T = -R * T / (V - b) ^ 2 + (2 * b + 2 * V) * (a * (1 + c) ^ 2 + 2 * a * (1 + c) * c * Tcm ^ -0.5 * T ^ 0.5 + a * c ^ 2 * Tcm ^ -1 * T) / (V ^ 2 + 2 * b * V - b ^ 2) ^ 2

            Cpm_ig = 0
            i = 0
            Do
                Cpm_ig += Vzmass(i) * VCpig(i)
                i += 1
            Loop Until i = n + 1

            Dim JT = -(T * dP_dT_V / dP_dV_T + V) / (Cpm_ig * MMm)

            JT_PR = JT*/
        }

        ///
        double Pvp_leekesler(double  T, double  Tc, double  Pc, double  w)
        {

            var Tr = T / Tc;

            var f0 = 5.92714 - 6.09648 / Tr - 1.28862 * Math.Log(Tr) + 0.169347 * Tr ^ 6
            var f1 = 15.2518 - 15.6875 / Tr - 13.4721 * Math.Log(Tr) + 0.43577 * Tr ^ 6

            var tmp = Pc * Math.Exp(f0 + w * f1);

            return  tmp;
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="mu_aq"></param>
        /// <param name="mu_ol"></param>
        /// <param name="xv_ol"></param>
        /// <returns></returns>
        double mu_ml(double mu_aq, double mu_ol, double xv_ol)
        {
            var mu_ml = 0.0;

            if( xv_ol >= 0.5)
            {
                mu_ml = mu_ol * Math.Exp(3.6 * (1 - xv_ol));
            }
            else if(xv_ol < 0.33 )
                mu_ml = (1 + 2.5 * xv_ol * (mu_ol + 0.4 * mu_aq) / (mu_ol + mu_aq)) * mu_aq;
            else
                mu_ml = xv_ol * mu_ol * Math.Exp(3.6 * (1 - xv_ol)) + (1 - xv_ol) * (1 + 2.5 * xv_ol * (mu_ol + 0.4 * mu_aq) / (mu_ol + mu_aq)) * mu_aq;

            return mu_ml;

        }
    }
}
